//
// Ephi - simulation of magnetic fields and particles
// Copyright (C) 2007 Indrek Mandre <indrek(at)mare.ee>
// For more information please see http://www.mare.ee/indrek/
//
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License along
// with this program; if not, write to the Free Software Foundation, Inc.,
// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
//

#ifndef __potential_hpp__
#define __potential_hpp__

#include <list>
#include <string>
#include <vector>

class PField;
class PLink;
class PMap;
class Statics;

class PPoint
{
public:
  PPoint (const vect3d& pos, const vect3d& ppos = vect3d(0, 0, 0), prec_t potential = 0);

  const vect3d& get_pos () const { return pos; }

private:
  typedef std::list<PLink*> links_t;
  links_t links;

  struct reflection_t
  {
    vect3d pos, pot_pos;
  };

  typedef std::vector<reflection_t> reflections_t;
  reflections_t reflections;

  size_t id;
  vect3d pos, pot_pos;
  prec_t potential;
  prec_t charge;
  prec_t rho;
  enum eflags { NO_CDIST = 1, IS_BOUNDARY = 2 };
  size_t flags;

  vect3d cur_efield;
  prec_t inc;

  void step() { charge += inc; inc = 0; }

  friend class PField;
  friend class PLink;
  friend class StaticPotentialField;
};

class PLink
{
public:
  PLink(PPoint *a, PPoint *b) : a(a), b(b), q(0)
  {
    ab = b->pos - a->pos;
    ab_normal = ab.normal();
    rr = ab*ab;
    len = prec_t_sqrt (rr);
    mid = 0.5 * (a->pos + b->pos);
    cc_normal = (0.05 * rr) * ab_normal;
    abr = (1 / (ab * ab)) * ab;
    rhodiff_by_llen = (b->rho - a->rho) / len;
    boundary = false;
    tmp = 0;
  }

private:
  PPoint *a, *b;
  vect3d ab;
  vect3d ab_normal;
  vect3d abr;
  vect3d cc_normal;
  prec_t rr;
  prec_t len;
  vect3d mid;
  prec_t rhodiff_by_llen;
  bool boundary;

  prec_t q;
  vect3d efield;
  PPoint *tmp;

  struct reflection_t
  {
    vect3d a, b;
    vect3d ab_normal;
  };
  typedef std::vector<reflection_t> reflections_t;
  reflections_t reflections;

  friend class PField;
  friend class PPoint;
  friend class StaticPotentialField;
};

class StaticPotentialField;

class PField
{
public:
  enum emethod_t { POINTSYSTEM = 1, CHARGEFLOW, LINESYSTEM };

  #define UNKNOWN_POTENTIAL 9e99

  PField (emethod_t method = POINTSYSTEM, prec_t resolution = 0.02);
  ~PField ();

  PPoint& addPoint (const vect3d& pos, const vect3d& ppos = vect3d(0, 0, 0), prec_t potential = 0);
  void addFeedPoint (PPoint& p, prec_t potential);

  void addSphere (const vect3d& p, prec_t radius, prec_t potential = 0, bool backawrds = false);
  void addTorus (const vect3d& p, const vect3d& normal, prec_t radius, prec_t wr, prec_t potential = 0);
  void addStrip (const vect3d& p0, const vect3d& p1, size_t n_segments = 0);
  void addT (const vect3d& p0, const vect3d& p1, const vect3d& p2, const vect3d& p3);
  void addSphereSegment (prec_t radius, const vect3d& n1, const vect3d& n2,
      const vect3d& n3, prec_t potential, bool outside = false, size_t level = 3);

  bool equalize();

  std::string getPovrayScript ();
  std::string getChargeConcentrationScript ();

  void getEfield(const vect3d& pos, vect3d& efield);
  prec_t getPotential(const vect3d& pos);

  void setPointCharge (const vect3d& pos, prec_t q);
  void setUniformCharge (const vect3d& pos, prec_t q, prec_t radius);

  void setStatics (Statics& statics) { this->statics = &statics; }

  PLink& connect (PPoint& a, PPoint& b);
  PLink* safe_connect (PPoint& a, PPoint& b, const std::vector<PLink*>& check);
  PLink& split (PLink& link, PPoint& mid);

  void splitSphericalTriangles (const vect3d& mid, prec_t radius,
      std::vector<PPoint*> spoints, std::vector<PLink*> slinks,
      prec_t potential, bool ouside = false, prec_t pratio = 0.0);

  void print();

  prec_t getAddedCharge () const { return added_charge; }
  prec_t getAveragePotential();
  void setPotential (prec_t pot);

  void save (const char *path);
  void save (const std::string& path) { save (path.c_str()); }
  bool load (const char *path);
  bool load (const std::string& path) { return load (path.c_str()); }

  void set_no_cdist (bool ncd) { no_cdist = ncd; }

  void getStatistics (size_t& pointcount, size_t& linkcount) { pointcount = points.size(); linkcount = links.size(); }

  prec_t getResolution() const { return resolution; }
  prec_t setResolution(prec_t res) { prec_t ores = resolution; resolution = res; return ores; }

  bool setTransforms (const std::vector<vect3d>& planes, const std::vector<transf3d>& transforms);

  bool fill (PMap& pmap);

  // splits the structures with double density
  void split ();

  emethod_t getMethod () const { return method; }

  bool empty() const { return points.empty(); }

private:
  emethod_t method;
  prec_t resolution;
  size_t pcounter;
  bool no_cdist;
  std::vector<vect3d> planes;
  std::vector<transf3d> transforms;

  typedef std::list<PPoint> points_t;
  points_t points;

  struct fpoint_t
  {
    fpoint_t() : cp(0), q(0) { }
    PPoint *point;
    prec_t potential;
    prec_t cp;
    prec_t q;
  };

  typedef std::list<fpoint_t> feeds_t;
  feeds_t feeds;

  prec_t added_charge;

  prec_t pc_q;
  vect3d pc_pos;
  prec_t pc_radius;
  prec_t pc_rr;
  prec_t pc_rho;

  Statics *statics;

  typedef std::list<PLink> links_t;
  links_t links;

  struct task : Task
  {
    PField *pf;
    links_t::iterator begin;
    links_t::iterator end;
    vect3d *cache;
    void execute () { pf->calc_efield (begin, end, cache); }
  };

  struct linesystem_task : Task
  {
    PField *pf;
    LinearSystem *ls;
    size_t part;
    size_t n;
    void execute () { pf->calc_linesystem (*ls, part, n); }
  };

  struct pointsystem_task : Task
  {
    PField *pf;
    LinearSystem *ls;
    size_t part;
    size_t n;
    void execute () { pf->calc_pointsystem (*ls, part, n); }
  };

  void reflect ();

  void calc_efield (const links_t::iterator& begin, const links_t::iterator& end, const vect3d* cache);
  void calc_linesystem (LinearSystem& ls, size_t part, size_t n);
  void calc_pointsystem (LinearSystem& ls, size_t part, size_t n);

  bool equalize_chargeflow();
  bool equalize_pointsystem();
  bool equalize_linesystem();

  PLink& unsafe_connect (PPoint& a, PPoint& b);

  vect3d * cache;

  friend class StaticPotentialField;
};

class StaticPotentialField : public StaticElement
{
public:
  StaticPotentialField (PField &pfield) : pfield(pfield) { }

  void addFields (const vect3d& pos, vect3d& bfield, vect3d& efield);
  void addPotential (const vect3d& pos, prec_t& potential);
  void draw(Screen& screen, const Screen::color& c = Screen::BLACK);
  std::string get_povray_script ();
  bool inContact (const vect3d& pos);

private:
  PField& pfield;
};

#endif // !__potential_hpp__

